library(tidyverse)
library(ggpubr)
library(raster)
library(ncdf4)
library(ozmaps)

###
# Functions needed for Thermal Sensitivity Perspective Code
###



###
# Find the intercept from GLM
findInt <- function(model, value) {
  function(x) {
    predict(model, data.frame(log10timeh = x), type = "response") - value
  }
}



###
# Calculate the accumulated damage based on time in minutes
# Function Source: Faber, Ørsted, Ehlers (2024) https://doi.org/10.1093/jxb/erae096 
Acc_function = function(x, y, intercept, slope) {
  output = data.frame(z = NA, cumsum = NA)
  last_row = nrow(data.frame(x = x, y = y))
  x = c(x, x[last_row])
  y = c(y, y[last_row])
  for (i in 1:length(x) - 1) {
    #print(paste0("round nr. ", i))
    acc = max(0, min(100, (100 * (x[i + 1] - x[i])) / (10^(slope * max(y[i + 1], y[i]) + intercept))))
    output[i, 1] = acc
  }
  output$cumsum = cumsum(output$z)
  z = output$cumsum
  z = append(z, 0)
  z[z > 100] = 100 # Set values exceeding 100% to 100%
  z = z[-length(z)]
  return(z)
}



###
# Repair function ArrFunc5 by Michael Kearney
ArrFunc5 <- function(x,TA,TAL,TAH,TL,TH,TREF,kdot_ref) {
  exp(TA * (1/TREF-1/(273.15+x)))/(1+exp(TAL * (1/(273.15+x)-1/TL)) + exp(TAH * (1/TH-1/(273.15+x)))) * kdot_ref
}



###
# Function to estimate thermal tolerance landscape from static assays
# General procedure
# Step 1: Calculate CTmax and z from TDT curve
# Step 2: Calculate average log10 time and Ta (mean x and y for interpolation purposes)
# Step 3: Interpolating survival probabilities to make them comparable across treatments 
# Step 4: Overlap all survival curves into a single one by shifting each curve to mean x and y employing z
# Step 5: Build expected survival curve with mean x and y pooling all data
# Step 6: Expand expected curve to multiple Ta (with 0.1ºC difference for predictive purposes)

tolerance.landscape <- function(ta,time){
  
  data <- data.frame(ta,time)
  data <- data[order(data$ta,data$time),]
  
  # Step 1: Calculate CTmax and z from TDT curve
  ta <- as.numeric(levels(as.factor(data$ta)))          
  model <- lm(log10(data$time) ~ data$ta); summary(model)
  ctmax <- -coef(model)[1]/coef(model)[2]
  z <- -1/coef(model)[2]
  
  # Step 2: Calculate average log10 time and Ta (mean x and y for interpolation purposes)
  time.mn <- mean(log10(data$time))
  ta.mn <- mean(data$ta)
  
  # Step 3: Interpolating survival probabilities to make them comparable across treatments 
  time.interpol <- matrix(NA,1001,length(ta))
  for(i in 1:length(ta)){   
    time <- c(0,sort(data$time[data$ta==ta[i]]))
    p <- seq(0,100,length.out = length(time))
    interp <- approx(p,time,n = 1001)
    time.interpol[,i] <- interp$y
  }         
  
  # Step 4: Overlap all survival curves into a single one by shifting each curve to mean x and y employing z
  # Step 5: Build expected survival curve with median survival time for each survival probability
  shift <- (10^((ta - ta.mn)/z))
  time.interpol.shift <- t(t(time.interpol)*shift)[-1,]
  surv.pred <- 10^apply(log10(time.interpol.shift),1,median)    
  
  # Step 6: Expand predicted survival curves to measured Ta (matrix m arranged from lower to higher ta)
  # Step 7: Obtain predicted values comparable to each empirical measurement
  m <- surv.pred*matrix ((10^((ta.mn - rep(ta, each = 1000))/z)), nrow = 1000)
  out <-0
  for(i in 1:length(ta)){
    time <- c(0,data$time[data$ta==ta[i]])
    p <- seq(0,100,length.out = length(time))
    out <- c(out,approx(seq(0,100,length.out = 1000),m[,i],xout=p[-1])$y)
  }
  data$time.pred <- out[-1]
  colnames(m) <- paste("time.at",ta,sep=".")
  m <- cbind(surv.prob=seq(1,0.001,-0.001),m)
  
  par(mfrow=c(1,2),mar=c(4.5,4,1,1),cex.axis=1.1)
  plot(-10,-10,las=1,xlab="Time (min)",ylab="Survival (%)",col="white",xaxs="i",yaxs="i",xlim=c(0,max(data$time)*1.05),ylim=c(0,105))
  for(i in 1:length(ta)){
    time <- c(0,sort(data$time[data$ta==ta[i]])); p <- seq(100,0,length.out = length(time))
    points(time,p,pch=21,bg="black",cex=0.5)
    time <- c(0,sort(data$time.pred[data$ta==ta[i]]))
    points(m[,i+1],100*m[,1],type="l",lty=2)}
  segments(max(data$time)*0.7,90,max(data$time)*0.8,90,lty=2)
  text(max(data$time)*0.82,90,"fitted",adj=c(0,0.5))
  plot(log10(data$time.pred),log10(data$time),pch=21,bg="black",cex=0.5,lwd=0.7,las=1,xlab="Fitted Log10 time",ylab="Measured Log10 time")
  abline(0,1,lty=2)
  rsq <- round(summary(lm(log10(data$time) ~ log10(data$time.pred)))$r.square,3)
  text(min(log10(data$time.pred)),max(log10(data$time)),substitute("r"^2*" = "*rsq),adj=c(0,1))
  list(ctmax = as.numeric(ctmax), z = as.numeric(z), ta.mn = ta.mn,  S = data.frame(surv=seq(0.999,0,-0.001),time=surv.pred),
       time.obs.pred=cbind(data$time,data$time.pred), rsq = rsq)
}                   



###
# Function dynamic.landscape1 by Michael Kearney and Enrico Rezende
dynamic.landscape1 <- function(ta,tolerance.landscape,kdot_ref){
  
  ArrFunc5 <- function(x,TA,TAL,TAH,TL,TH,TREF,kdot_ref){
    exp(TA*(1/TREF-1/(273.15+x)))/(1+exp(TAL*(1/(273.15+x)-1/TL))+exp(TAH*(1/TH-1/(273.15+x))))*kdot_ref
  }
  
  TA <- 1000
  TL <- 8.5 + 273.15
  TH <- 25.5 + 273.15
  TAL <- 100000
  TAH <- 100000
  TREF <- 20 + 273.15
  
  Tbs <- seq(0, 50)
  TempCorr <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, 0.001)
  
  surv <- tolerance.landscape$S[,2] + seq(0,0.0000001,length.out=1000)
  ta.mn <- tolerance.landscape$ta.mn
  z <- tolerance.landscape$z
  shift <- 10^((ta.mn - ta)/z)  
  time.rel <- 0
  alive <- 100
  survival <- matrix(NA, length(ta))
  repair <- ArrFunc5(ta, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
  for(i in 1:length(ta)){           
    if(!is.na(alive)) {
      alive <- try(approx(c(0,shift[i]*surv),seq(100,0,length.out = length(c(0,surv))),
                          xout = time.rel[i] + 1)$y,silent=TRUE)    
      alive <- alive + ArrFunc5(ta[i], TA, TAL, TAH, TL, TH, TREF, kdot_ref)
      alive <- ifelse(alive > 100,100,alive)
      survival[i] <- alive
      time.rel <- try(c(time.rel,approx(seq(100,0,length.out = length(c(0,surv))),c(0,shift[i + 1]*surv),
                                        xout = alive)$y),silent=TRUE)}
    else{
      alive <- 100
      survival[i] <- 0}}                
  out <- data.frame(cbind(ta=ta[1:(length(survival)-1)],time=(1:length(ta))[1:(length(survival)-1)],
                          alive=survival[1:(length(survival)-1)]))
  par(mar=c(4,4,1,1),mfrow=c(1,2))
  plot(1:length(ta),ta,type="l",xlim=c(0,length(ta)),ylim=c(min(ta),max(ta)),col="black",lwd=1.5,las=1,
       xlab = "Time (min)", ylab = "Temperature (ºC)")          
  plot(out$time,out$alive,type="l",xlim=c(0,length(ta)),ylim=c(0,100),col="black",lwd=1.5,las=1,
       xlab = "Time (min)", ylab = "Survival (%)")
  list(time = out$time,ta = out$ta, survival = out$alive)}


# Add in a dynamic repair function within by Enrico Rezende, Michael Kearney, Pieter Arnold
dynamic.landscape2 <-
  function(ta,tolerance.landscape,kdot_ref){
  
  ArrFunc5 <- function(x,TA,TAL,TAH,TL,TH,TREF,kdot_ref){
    exp(TA*(1/TREF-1/(273.15+x)))/(1+exp(TAL*(1/(273.15+x)-1/TL))+exp(TAH*(1/TH-1/(273.15+x))))*kdot_ref
  }
  
  TA <- 1000
  TL <- 8.5 + 273.15
  TH <- 25.5 + 273.15
  TAL <- 100000
  TAH <- 100000
  TREF <- 20 + 273.15
  
  Tbs <- seq(0, 50)
  TempCorr <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, 0.001)
  
  surv <- tolerance.landscape$S[,2] + seq(0,0.0000001,length.out=1000)
  ta.mn <- tolerance.landscape$ta.mn
  z <- tolerance.landscape$z
  shift <- 10^((ta.mn - ta)/z)  
  time.rel <- 0
  alive <- 100
  survival <- matrix(NA, length(ta))
  kdot <- matrix(NA, length(ta))
  repair <- ArrFunc5(ta, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
  for(i in 1:length(ta)){           
    if(!is.na(alive)) {
      alive <- try(approx(c(0,shift[i]*surv),seq(100,0,length.out = length(c(0,surv))),
                          xout = time.rel[i] + 1)$y,silent=TRUE)    
      # modify kdot_ref based on alive. If alive > 99, kdot_ref stays the same. If alive < 99, kdot_ref reduces as a function of alive
      kdot_ref1 <- ifelse(alive < 99, kdot_ref*(kdot_ref^(3/(alive))), kdot_ref)
      alive <- alive + ArrFunc5(ta[i], TA, TAL, TAH, TL, TH, TREF, kdot_ref1)
      alive <- ifelse(alive > 100,100,alive)
      survival[i] <- alive
      kdot[i] <- kdot_ref1
      time.rel <- try(c(time.rel,approx(seq(100,0,length.out = length(c(0,surv))),c(0,shift[i + 1]*surv),
                                        xout = alive)$y),silent=TRUE)}
    else{
      alive <- 100
      survival[i] <- 0
      kdot[i] <- kdot_ref }}                
  out <- data.frame(cbind(ta=ta[1:(length(survival)-1)],time=(1:length(ta))[1:(length(survival)-1)],
                          alive=survival[1:(length(survival)-1)],kdot=kdot[1:(length(survival)-1)]))
  par(mar=c(4,4,1,1),mfrow=c(1,2))
  plot(1:length(ta),ta,type="l",xlim=c(0,length(ta)),ylim=c(min(ta),max(ta)),col="black",lwd=1.5,las=1,
       xlab = "Time (min)", ylab = "Temperature (ºC)")          
  plot(out$time,out$alive,type="l",xlim=c(0,length(ta)),ylim=c(0,100),col="black",lwd=1.5,las=1,
       xlab = "Time (min)", ylab = "Survival (%)")
  list(time = out$time,ta = out$ta, survival = out$alive, kdot = out$kdot)}




###
# Function to extract data from netcdf files by Michael Kearney
extract.nc <- function(file, loc){
  nc <- nc_open(file)
  lon <- ncvar_get(nc, "longitude") # extract longitude values
  lat <- ncvar_get(nc, "latitude") # extract latitude values
  # find closest pixel to site requested
  flat <- match(abs(lat-loc[2])<.3,1)
  latindex <- which(flat %in% 1)
  flon <- match(abs(lon-loc[1])<.3,1)
  lonindex <- which(flon %in% 1)
  start <- c(1,latindex,lonindex) # start location in file for extraction
  count <- c(-1, 1, 1) # extract all data (-1)
  data <- as.numeric(ncvar_get(nc, varid = attributes(nc$var)$names, start = start, count)) # get the data
  # extract dates
  dates <- as.POSIXct(ncvar_get(nc, "time"), tz = "Etc/GMT+10", origin = "1970-01-01")
  nc_close(nc)
  result <- cbind(dates, as.data.frame(data)) # final output
}


#source("Thermal sensitivity functions.R")

cbp2 <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
          "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
weeddata <- read.csv("weedseeddata_all.csv") 
# Data extracted from Dahlquist et al. (2007) Figures using metaDigitise
# https://doi.org/10.1614/WS-04-178.1 

weeddata$mortality_prop <- weeddata$mortality_percent/100
weeddata$log10timemin <- log10(weeddata$time_min)
weeddata$log10timeh <- log10(weeddata$time_h)
weeddata_sp1 <- subset(weeddata, species=="London rocket")

# Load in extracted dataset that contains LT80 values
weedLT80 <- read.csv("Weed_LT80.csv")
weedLT80$log10timemin <- log10(weedLT80$time_min)
weedLT80$log10timeh <- log10(weedLT80$time_h)
model50 <- glm(mortality_prop ~ log10timeh, 
             data = weeddata_sp1[weeddata_sp1$temp_treatment=="50",], family = quasibinomial)
model46 <- glm(mortality_prop ~ log10timeh, 
             data = weeddata_sp1[weeddata_sp1$temp_treatment=="46",], family = quasibinomial)
model42 <- glm(mortality_prop ~ log10timeh, 
             data = weeddata_sp1[weeddata_sp1$temp_treatment=="42",], family = quasibinomial)
# Fit overall and species by species TDT models
overallLT80_mod <- lm(temp_treatment ~ log10timeh, data = weedLT80)

LT80mod <- lm(temp_treatment ~ log10timeh, data = weedLT80)
LT80mod_sp3 <- lm(temp_treatment ~ log10timeh, data = subset(weedLT80, weedLT80$Species2=="Black nightshade"))
LT80mod_sp5 <- lm(temp_treatment ~ log10timeh, data = subset(weedLT80, weedLT80$Species2=="London rocket"))
LT80mod_sp6 <- lm(temp_treatment ~ log10timeh, data = subset(weedLT80, weedLT80$Species2=="Tumble pigweed"))

# Fit the alternative model form: time ~ temp to apply accumulated damage model
LT80mod_Orsted <- lm(log10timeh ~ temp_treatment, data = weedLT80)
LT80mod_Orsted1 <- lm(log10timeh ~ temp_treatment, data = subset(weedLT80, weedLT80$Species2=="Black nightshade"))
LT80mod_Orsted2 <- lm(log10timeh ~ temp_treatment, data = subset(weedLT80, weedLT80$Species2=="London rocket"))
LT80mod_Orsted3 <- lm(log10timeh ~ temp_treatment, data = subset(weedLT80, weedLT80$Species2=="Tumble pigweed"))

intercept <- as.numeric(LT80mod_Orsted1$coefficients[1])
slope <- as.numeric(LT80mod_Orsted1$coefficients[2])

colfunc <- colorRampPalette(c("palegreen","tan4"))
facet.labs <- c("Solanum nigrum", "Sisymbrium irio", "Amaranthus albus")
names(facet.labs) <- c("Acc_damage_sp1", "Acc_damage_sp2", "Acc_damage_sp3")

# Simulate 50ËšC temperature ramp and hold treatment with thermal variability
#st1 <- c(seq(30,39,0.5),rep(40,20),seq(41,45,0.5),rep(46,20),seq(47,51,0.5),rep(52,100))
st1 <- c(seq(20,50,2.5),rep(50,100))
specific_temperatures1 <- st1 + rnorm(st1,0,1)
Temperature <- specific_temperatures1
Time <- specific_temperatures1
data <- data.frame(Time,Temperature)
data <- data %>% filter(row_number() %% 2 != 1)
data$Time <- seq(1, nrow(data))
#data$Time <- as.numeric(data$Time)/2
data$Temperature <- as.numeric(data$Temperature)
data <- data[1:100,]

# Calculate accumulated damage
data$Acc_damage_mean <- Acc_function(data$Time,data$Temperature, 
                                as.numeric(LT80mod_Orsted$coefficients[1]), # intercept
                                as.numeric(LT80mod_Orsted$coefficients[2])) # slope 
data$Acc_damage_mean <- as.numeric(data$Acc_damage_mean)
data$Acc_damage_sp1 <- Acc_function(data$Time,data$Temperature, 
                                as.numeric(LT80mod_Orsted1$coefficients[1]), # intercept
                                as.numeric(LT80mod_Orsted1$coefficients[2])) # slope 
data$Acc_damage_sp2 <- Acc_function(data$Time,data$Temperature, 
                                as.numeric(LT80mod_Orsted2$coefficients[1]), # intercept
                                as.numeric(LT80mod_Orsted2$coefficients[2])) # slope 
data$Acc_damage_sp3 <- Acc_function(data$Time,data$Temperature, 
                                as.numeric(LT80mod_Orsted3$coefficients[1]), # intercept
                                as.numeric(LT80mod_Orsted3$coefficients[2])) # slope 

datalong <- data %>% pivot_longer(cols = c(Acc_damage_sp1, Acc_damage_sp2, Acc_damage_sp3), 
                                  names_to = "Acc_damage")
datalong <- na.omit(datalong)
preddata <- matrix(data = NA, 51, 3)
preddata[,2] <- seq(0,5,0.1)
preddata <- as.data.frame(preddata)
colnames(preddata) <- c("test", "log10timeh","out")

# CTmax1h temperature (at log10(0))
predict(LT80mod, as.data.frame(preddata))[1]
##        1 
## 57.13306
LT80mod$coefficients[1] # CTmax1h equivalent (just checking the predict works)
## (Intercept) 
##    57.13306
LT80mod_sp3$coefficients[1]
## (Intercept) 
##    60.09522
LT80mod_sp5$coefficients[1]
## (Intercept) 
##    53.47222
LT80mod_sp6$coefficients[1]
## (Intercept) 
##    59.32692
# Find time for LT80 in log10 hours
uniroot(findInt(model50, 0.8), range(weeddata_sp1$log10timeh))$root
## [1] 0.60383
uniroot(findInt(model46, 0.8), range(weeddata_sp1$log10timeh))$root
## [1] 1.305468
uniroot(findInt(model42, 0.8), range(weeddata_sp1$log10timeh))$root
## [1] 1.929779
# Find time for LT80 in hours
10^(uniroot(findInt(model50, 0.8), range(weeddata_sp1$log10timeh))$root)
## [1] 4.016336
10^(uniroot(findInt(model46, 0.8), range(weeddata_sp1$log10timeh))$root)
## [1] 20.2054
10^(uniroot(findInt(model42, 0.8), range(weeddata_sp1$log10timeh))$root)
## [1] 85.07045
# Hours at treatment of 50ËšC to reach LT80
10^(uniroot(findInt(LT80mod, 50), range(weedLT80$log10timeh), extendInt = "yes")$root)     # Average
## [1] 19.68409
10^(uniroot(findInt(LT80mod_sp3, 50), range(weedLT80$log10timeh), extendInt = "yes")$root) # Black nightshade
## [1] 35.07929
10^(uniroot(findInt(LT80mod_sp5, 50), range(weedLT80$log10timeh), extendInt = "yes")$root) # London rocket
## [1] 3.719188
10^(uniroot(findInt(LT80mod_sp6, 50), range(weedLT80$log10timeh), extendInt = "yes")$root) # Tumble pigweed
## [1] 61.95814
# LT80 for london rocket seeds

fig1A <- 
ggplot(data = weeddata_sp1, aes(y = mortality_prop, x = log10timeh)) +
  geom_point(aes(colour = as.factor(temp_treatment)), alpha = 0.4) +
  geom_smooth(aes(colour = as.factor(temp_treatment), fill = as.factor(temp_treatment)), 
              method = "glm", formula = y ~ x,
    method.args = list(family = "quasibinomial"), se = T, fullrange = T) +
  geom_hline(yintercept = 0.8, lty = 2) +
  geom_segment(x = 1.929779, xend = 1.929779, y = 0.8, yend = -Inf, lty = 2) + # 42
  geom_segment(x = 1.305468, xend = 1.305468, y = 0.8, yend = -Inf, lty = 2) + # 46
  geom_segment(x = 0.60383, xend = 0.60383, y = 0.8, yend = -Inf, lty = 2) + # 50
  scale_colour_manual(values = c("goldenrod1", "darkorange2", "red3"), "Temperature (°C)") +
  scale_fill_manual(values = c("goldenrod1", "darkorange2", "red3"), "Temperature (°C)") +
  annotate("text", x = 0.4, y = 1.1, label = bquote(italic(Sisymbrium~irio))) +
  # annotate("text", x = 0.4, y = 0.58, label = "50°C") +
  # annotate("text", x = 1.1, y = 0.58, label = "46°C") +
  # annotate("text", x = 1.7, y = 0.58, label = "42°C") +
  annotate("text", x = 0, y = 0.75, label = bquote(italic(LT)[80])) +
  annotate("text", x = 2.5, y = 0.75, label = bquote(italic(LT)[80])) +
  labs(y = "Seed mortality (proportion)", x = bquote(log[10]~Time~(h))) +
  scale_x_continuous(limits = c(-0.4,3)) +
  scale_y_continuous(limits = c(0,1.2), breaks = c(0,0.2,0.4,0.6,0.8,1)) +
  theme_classic() + theme(legend.position = "bottom", legend.box = "vertical")

#fig1A

# TDT models fit to 3 weed species and overall, CTmax' and z

fig1B <-
ggplot(data = weedLT80, aes(y = temp_treatment, x = log10(time_h), 
                                  group = Species, fill = Species, shape = Species)) +
  geom_abline(intercept = LT80mod_sp3$coefficients[1], slope = LT80mod_sp3$coefficients[2], 
              colour = cbp2[6], lty = 1, linewidth = 1) + # black_nightshade, solanum nigrum
  geom_abline(intercept = LT80mod_sp5$coefficients[1], slope = LT80mod_sp5$coefficients[2], 
              colour = cbp2[4], lty = 1, linewidth = 1) + # london_rocket, sisymbrium irio
  geom_abline(intercept = LT80mod_sp6$coefficients[1], slope = LT80mod_sp6$coefficients[2], 
              colour = cbp2[8], lty = 1, linewidth = 1) + # tumble_pigweed, amaranthus albus
  #geom_abline(intercept = overallLT80_mod$coefficients[1], slope = overallLT80_mod$coefficients[2], 
  #            colour = cbp2[1], lty = 1, linewidth = 1) + # overall
  geom_point(size = 3) +
  geom_vline(xintercept = 0, lty = 2) +
  scale_fill_manual(values = c(cbp2[6], cbp2[4], cbp2[8], cbp2[1]), name = "",
                     labels = c(bquote(italic(Solanum~nigrum)), 
                                bquote(italic(Sisymbrium~irio)), 
                                bquote(italic(Amaranthus~albus)))) +
  scale_shape_manual(values = c(22,23,24), name = "",
                     labels = c(bquote(italic(Solanum~nigrum)), 
                                bquote(italic(Sisymbrium~irio)), 
                                bquote(italic(Amaranthus~albus)))) +
  annotate("text", x = 0.2, y = 60, label = bquote(italic(CT[max1h])~"="~intercept), hjust = 0) +
  annotate("text", x = 1.6, y = 54, label = bquote(italic(z)~"="~slope)) +
  annotate("text", x = 0.05, y = 63, label = bquote(italic(LT)[80]~TDT~curve), hjust = 0) +
  annotate("text", x = 0.05, y = 42.5, label = bquote(italic(CT[max1h])~"="~"60.1°C"), colour = cbp2[6], hjust = 0) +
  annotate("text", x = 0.05, y = 41, label = bquote(italic(CT[max1h])~"="~"53.5°C"), colour = cbp2[4], hjust = 0) +
  annotate("text", x = 0.05, y = 39.5, label = bquote(italic(CT[max1h])~"="~"59.3°C"), colour = cbp2[8], hjust = 0) +
  annotate("text", x = 0.05, y = 38, label = bquote(italic(z)==-6.534), colour = cbp2[6], hjust = 0) +
  annotate("text", x = 0.05, y = 36.5, label = bquote(italic(z)==-6.087), colour = cbp2[4], hjust = 0) +
  annotate("text", x = 0.05, y = 35, label = bquote(italic(z)==-5.205), colour = cbp2[8], hjust = 0) +
  scale_x_continuous(limits = c(-0.2,3)) +
  scale_y_continuous(limits = c(35,65), breaks = c(35,40,45,50,55,60,65)) +
  labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
  theme_classic() + theme(legend.position = "bottom")

#fig1B

# Damage accumulation model applied to 3 species over temperature-time treatment simulation

fig1C <- 
ggplot(datalong, aes(x = Time, y = Temperature)) +
  facet_wrap(~Acc_damage, nrow = 3, labeller = labeller(Acc_damage = facet.labs)) +
  geom_line() +
  geom_point(data = datalong, aes(x = Time, y = Temperature, fill = value),
             shape = 21, size = 1.5, alpha = 1) +
  scale_fill_gradientn("Damage (%)", colours = colfunc(21), limits = c(0,100)) +
  labs(x = "Time (h)", y = bquote(Temperature~(degree*C))) +
  scale_x_continuous(limits = c(0,60), breaks = c(0,10,20,30,40,50)) +
  ylim(20,55) +  
  theme_classic() + 
  theme(legend.position = "right", strip.background = element_blank(),
        strip.text = element_text(face = "italic"))
#fig1C
  
# Damage accumulation model applied to 3 species over temperature-time treatment simulation

fig1D <- 
  ggplot(data, aes(x = Time, y = Acc_damage_mean)) +
  geom_line(data = data, aes(x = Time, y = Acc_damage_sp1, colour = cbp2[6]), linewidth = 1) +
  geom_line(data = data, aes(x = Time, y = Acc_damage_sp2, colour = cbp2[4]), linewidth = 1) +
  geom_line(data = data, aes(x = Time, y = Acc_damage_sp3, colour = cbp2[8]), linewidth = 1) +
  #geom_point(data = data, aes(x = Time, y = Acc_damage_sp1, fill = cbp2[6]), shape = 22, size = 2, alpha = 0.5) +
  #geom_point(data = data, aes(x = Time, y = Acc_damage_sp2, fill = cbp2[4]), shape = 23, size = 2, alpha = 0.5) +
  #geom_point(data = data, aes(x = Time, y = Acc_damage_sp3, fill = cbp2[8]), shape = 24, size = 2, alpha = 0.5) +
  scale_fill_manual(values = c(cbp2[6], cbp2[4], cbp2[8]), name = "",
                     labels = c(bquote(italic(Solanum~nigrum)), 
                                bquote(italic(Sisymbrium~irio)), 
                                bquote(italic(Amaranthus~albus)))) +
  scale_colour_manual(values = c(cbp2[6], cbp2[4], cbp2[8]), name = "",
                     labels = c(bquote(italic(Solanum~nigrum)), 
                                bquote(italic(Sisymbrium~irio)), 
                                bquote(italic(Amaranthus~albus)))) +
  scale_y_continuous(limits = c(0,110), breaks = c(0,25,50,75,100)) +
  scale_x_continuous(limits = c(0,60), breaks = c(0,10,20,30,40,50)) +
  geom_hline(yintercept = 100, linetype = "dashed") +
  labs(x = "Time (h)", y = expression(paste("Accumulated damage"," (%)"))) +
  theme_classic() + theme(legend.position = "none")

#fig1D

ggarrange(fig1A, fig1B, fig1C, fig1D,
          nrow = 2, ncol = 2, align = "v", labels = c("(a)", "(b)", "(c)", "(d)"))

#ggsave("Fig1_R1.pdf", height = 8.5, width = 10)
set.seed(21)

# Set Arrhenius model and repair rate parameters        
TA <- 14065 # Arrhenius temperature, K
TL <- 10.5 + 273.15 # Arrhenius temperature lower threshold, K
TH <- 28.5 + 273.15 # Arrhenius temperature upper threshold, K
TAL <- 50000 # Arrhenius temperature lower, K
TAH <- 100000 # Arrhenius temperature upper, K
TREF <- 20 + 273.15 # reference temperature, K
Tbs <- seq(0, 50) # sequence of temperatures over which to model, deg C
kdot_ref <- 0.02 # repair rate % per minute at ref temp 20 deg C

# Load ectotherm model data and subset to example four week period
par(mfrow = c(1,1))
ta <- read.csv("ectotherm.csv")[3006 + 1:(24*28),] # subset ectotherm - derived from NicheMapR
ta <- ta$TC + 0.3 # Body temperature
plot(ta, type = "l")

time.min <- spline(ta, n = length(ta-1)*60)$x   
ta.min <- spline(ta,n = length(ta-1)*60)$y  
ta <- ta.min
ind <- 1000
ctmax <- 43.8 # CTmax1h
z <-  -2.4 # z (slope)

static <- rep(c(round(ctmax,2),round(ctmax+z,2), round(ctmax+z*2,2)), each = ind)
time <- abs(c(rnorm(ind,1,1/4), rnorm(ind,10,10/4), rnorm(ind,100,100/4)))  
tl <- tolerance.landscape(static, time)

# Fit modified dynamic.landscape function that includes the Arrhenius repair function

# Set repair coefficients
kdot0 <- 0 # none
kdot1 <- 0.007 # low
kdot2 <- 0.0111 # moderate
kdot3 <- 0.02 # high

# No repair
dl0 <- dynamic.landscape2(ta,tl,kdot0)

# Three levels of repair rates
dl1 <- dynamic.landscape2(ta,tl,kdot1)

dl2 <- dynamic.landscape2(ta,tl,kdot2)

dl3 <- dynamic.landscape2(ta,tl,kdot3)

# Temperature-time model sequence
time2 <- log10(seq(1,330,1))
ht <- 43.8 + (-2.4*time2)
sim <- data.frame(cbind(time2,ht))
Tbs <- seq(5, 45, 0.2)

# Generate Arrhenius functions for repair rates
TempCorr0 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot0)
TempCorr1 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot1)
TempCorr2 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot2)
TempCorr3 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot3)

# Combine dataset
damage <- c(seq(0.000001,0.005,0.00005),(seq(0.1,3,0.01)^2.3))/15
damage <- damage[1:201]
range(damage)
## [1] 6.666667e-08 8.300646e-02
TempCorrDat <- data.frame(cbind(Tbs, TempCorr0, TempCorr1, TempCorr2, TempCorr3, damage))
TempCorrDatL <- TempCorrDat %>% pivot_longer(cols = c(TempCorr0, TempCorr1, TempCorr2, TempCorr3), 
                                             names_to = c("repair"), values_to = "repair_value")
TempCorrDatL <- subset(TempCorrDatL, Tbs > 6)
TempCorrDatLrepair <- subset(TempCorrDatL, repair!="TempCorr0") 

dldat <- data.frame(cbind(data.frame(dl0),data.frame(dl1),data.frame(dl2),data.frame(dl3)))
dldatL <- dldat %>% pivot_longer(cols = c(survival, survival.1, survival.2, survival.3), names_to = c("repair"),
                                 values_to = "repair_value")

dldatKL <- dldat %>% pivot_longer(cols = c(kdot, kdot.1, kdot.2, kdot.3), names_to = c("kdot"), 
                                 values_to = "kdot_value")

dldat_comb <- data.frame(cbind(dldatL$repair, dldatL$repair_value, dldatKL$kdot, dldatKL$kdot_value))
colnames(dldat_comb) <- c("repair", "repair_value", "kdot", "kdot_value")
dldat_comb$repair_value <- as.numeric(dldat_comb$repair_value)
dldat_comb$kdot_value <- as.numeric(dldat_comb$kdot_value)

func_seq <- seq(0,100,1)
repair0_func <- kdot0*(kdot0^(3/(func_seq)))
repair1_func <- kdot1*(kdot1^(3/(func_seq)))
repair2_func <- kdot2*(kdot2^(3/(func_seq)))
repair3_func <- kdot3*(kdot3^(3/(func_seq)))
rate_decay <- data.frame(func_seq, repair0_func, repair1_func, repair2_func, repair3_func)
rate_decayL <- rate_decay %>% pivot_longer(cols = c(repair0_func, repair1_func, repair2_func, repair3_func), 
                                           names_to = "repair", values_to = "repair_value")

# Plot sensitivity temp~time and display CTmax and z
base_tls_plot <- 
ggplot(data = sim, aes(y = ht, x = time2)) +
  geom_smooth(method = "lm", formula = y~x, colour = "black", fullrange = F) + 
  ylim(35,45) + xlim(0,2.6) +
  labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
  annotate("text", x = 0.3, y = 38, label = bquote(italic(CT[max1h])=="43.8°C"), hjust = 0) +
  annotate("text", x = 0.3, y = 36, label = bquote(italic(z)=="-2.4"), hjust = 0) +
  theme_classic()

# Damage rate
damage_plot <- 
ggplot(data.frame(TempCorrDatL), aes(x = Tbs, y = damage)) +
  geom_path() +
  scale_x_continuous(limits = c(5,45)) +
  scale_y_continuous(limits = c(0,0.06)) +
  labs(y = bquote(Damage~rate~("%"~min^-1)), x = bquote(Temperature~(degree*C))) +
  theme_classic()

# Repair rate
repair_plot <- 
ggplot(data.frame(TempCorrDatL), aes(x = Tbs, y = repair_value, colour = repair)) +
  geom_path(aes(group = rev(repair))) +
  scale_x_continuous(limits = c(5,45)) +
  scale_y_continuous(limits = c(0,0.06)) +
  scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]), 
                      labels = c("None", "Low", "Moderate", "High")) +
  labs(y = bquote(Repair~rate~("%"~min^-1)), x = bquote(Temperature~(degree*C))) +
  theme_classic()

# Ratio of damage:repair
ratio_plot <- 
ggplot(data.frame(TempCorrDatLrepair), aes(x = Tbs, y = (damage/repair_value), colour = repair)) +
  geom_path(aes(group = rev(repair))) +
  geom_hline(yintercept = 1, linetype = 2) +
  scale_x_continuous(limits = c(5,45)) +
  scale_y_log10(limits = c(0.005,10000), breaks = c(0.01, 0.1, 1, 10, 100, 1000), 
                labels = c("0.01", "0.1", "1", "10", "100", "1000")) +
  scale_colour_manual(name = "Repair rate", values = c(cbp2[2:4]), 
                      labels = c("Low", "Moderate", "High")) +
  annotate("text", x = 32, y = 0.05, label = "Repair \noutweighs \ndamage", hjust = 0) +
  annotate("text", x = 12, y = 1000, label = "Damage \noutweighs \nrepair", hjust = 0) +
  labs(y = bquote(Damage:repair~ratio), x = bquote(Temperature~(degree*C))) +
  theme_classic() + theme(legend.position = "none")

# Net damage rate
net_plot <- 
ggplot(data.frame(TempCorrDatL), aes(x = Tbs, y = damage-repair_value, colour = repair)) +
  geom_path(aes(group = rev(repair))) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_x_continuous(limits = c(5,45)) +
  scale_y_continuous(limits = c(-0.06,0.1)) +
  scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]), 
                      labels = c("None", "Low", "Moderate", "High")) +
  annotate("text", x = 32, y = -0.03, label = "Repair \noutweighs \ndamage", hjust = 0) +
  annotate("text", x = 12, y = 0.08, label = "Damage \noutweighs \nrepair", hjust = 0) +
  labs(y = bquote(Net~damage~rate~("%"~min^-1)), x = bquote(Temperature~(degree*C))) +
  theme_classic()

# Plot decline in repair over time due to damage accumulation 
decline_plot <- 
ggplot(dldatKL, aes(y = kdot_value, x = time, colour = kdot)) +
  geom_vline(xintercept = 6350, colour = "red", lty = 2) + 
  geom_vline(xintercept = 30870, colour = "red", lty = 2) + 
  geom_vline(xintercept = 36550, colour = "red", lty = 2) + 
  geom_path(aes(group = kdot)) +
  #scale_y_continuous(limits = c(0,100)) +
  scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]), 
                      labels = c("None", "Low", "Moderate", "High")) +
  labs(x = "Time (days)", y = bquote(Repair~rate~coefficient~(italic(dot(k))))) +
  theme_classic() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

# Plot four weeks of body temperature
tb_plot <- 
ggplot(data.frame(dl0), aes(y = ta, x = time)) +
  geom_hline(yintercept = 28.5, colour = "black", lty = 2) + 
  geom_vline(xintercept = 6350, colour = "red", lty = 2) + 
  geom_vline(xintercept = 30870, colour = "red", lty = 2) + 
  geom_vline(xintercept = 36550, colour = "red", lty = 2) + 
  geom_line() +
  labs(x = "Time (days)", y = bquote(Temperature~(degree*C))) +
  scale_y_continuous(limits = c(3,40)) +
  annotate("text", x = 12000, y = 38, label = bquote(italic(T)[c]=="28.5°C"), hjust = 0) +
  theme_classic() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

# Plot physiological function over the four weeks incorporating repair 
func_plot <- 
ggplot(dldatL, aes(y = repair_value, x = time, colour = repair)) +
  geom_vline(xintercept = 6350, colour = "red", lty = 2) + 
  geom_vline(xintercept = 30870, colour = "red", lty = 2) + 
  geom_vline(xintercept = 36550, colour = "red", lty = 2) + 
  geom_path(aes(group = rev(repair))) +
  scale_y_continuous(limits = c(0,100)) +
  scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]), 
                      labels = c("None", "Low", "Moderate", "High")) +
  labs(x = "Time (days)", y = "Physiological function (%)") +
  theme_classic() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

# Plot decline as a function of damage accumulation
decline_func_plot <- 
ggplot(rate_decayL, aes(y = repair_value, x = func_seq, colour = repair)) +
  geom_path(aes(group = rev(repair))) +
  #scale_y_continuous(limits = c(0,100)) +
  scale_x_reverse() +
  scale_colour_manual(name = "Repair rate", values = c(cbp2[1:4]), 
                      labels = c("None", "Low", "Moderate", "High")) +
  labs(x = "Physiological function (%)", y = bquote(Repair~rate~coefficient~(italic(dot(k))))) +
  theme_classic()

ggarrange(base_tls_plot, damage_plot, repair_plot, 
          ratio_plot, net_plot, decline_func_plot, 
          tb_plot, func_plot, decline_plot,
          labels = c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)", "(g)", "(h)", "(i)"), 
          align = "hv", ncol = 3, nrow = 3,
          common.legend = TRUE, legend = "bottom")

#ggsave("Fig2_R1.pdf", width = 10, height = 10)
# Load in climate data and crown dieback data
env <- read.csv("Summer_Climate_Penrith.csv", header = TRUE)
env$date <- dmy(as.character(env$Date)) 
env.noNA <- env[!is.na(env$Rainfall_mm),] # remove missing values from one column
scaleFactor <- max(env.noNA$Max_T_C) / max(env.noNA$Rainfall_mm)
die <- read.csv("Crown_Dieback.csv", header = TRUE)
die$Date <- dmy(as.character(die$Date))

fig3a <- ggplot(env, aes(x = date)) + 
  geom_rect(aes(NULL, NULL, xmin = as.Date("2019-12-5"), xmax = as.Date("2020-1-15")), ymin = 45, ymax = 55, 
            fill = alpha("burlywood", 0.1)) + 
  geom_rect(aes(NULL, NULL, xmin = as.Date("2020-1-16"), xmax = as.Date("2020-2-10")), ymin = 45, ymax = 55, 
            fill = alpha("lightblue", 0.1)) + 
  geom_rect(aes(NULL, NULL, xmin = as.Date("2020-2-11"), xmax = as.Date("2020-4-1")), ymin = 45, ymax = 55, 
            fill = alpha("seagreen3", 0.05)) +
  geom_line(aes(y = Max_T_C), colour = "black", linewidth = 0.5) +
  geom_bar(aes(y = Rainfall_mm * scaleFactor), fill = alpha("blue", 0.4), 
           stat = "identity", width = 0.8, size = 0.25) +
  coord_cartesian(xlim=c(as.Date("2019-12-15"), as.Date("2020-03-26")), ylim = c(0, 56)) +
  scale_x_date(date_labels = "%d-%b-%y") +
  scale_y_continuous(expand = c(0, 0), 
                     limits = c(0, 56),
                     sec.axis = sec_axis(~./scaleFactor, name="Rainfall (mm)")) +
  labs(x = "Date", y = expression("Daily Maximum"~italic("T")[air]~~(degree*C))) +
  theme(axis.title.x = element_blank(), axis.text.x = element_blank()) +
  annotate("text", x = as.Date("2019-12-29"), y = 51, 
           label = "Too dry to repair \n heat stress \n damage", 
           size = 3, colour = "orangered") +
  annotate("text", x = as.Date("2020-1-28"), y = 51, 
           label = "Unfavourable \n conditions \n for repair", 
           size = 3, colour = "blue3") + 
  annotate("text", x = as.Date("2020-3-5"), y = 51, 
           label = "Conditions allow repair \n of damage and canopy \n dieback regeneration", 
           size = 3, colour = "darkgreen") +
  theme_classic()

fig3a

#ggsave("Fig3a_R1.pdf", height = 5, width = 5)

fig3b <- 
ggplot(die, aes(x = Date, y = Survival_percent, colour = Species_Code)) + 
  geom_point(size = 2) +
  geom_line(aes(group = Species_Code)) +
  coord_cartesian(xlim = c(as.Date("2019-12-15"), as.Date("2020-03-26")), ylim = c(0, 100)) +
  scale_x_date(date_labels = "%d-%b-%y") +
  scale_colour_manual(values = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3")) +
  labs(x="Date", y="Tree crown survival (%)") + 
  annotate("text", x = as.Date("2020-1-12"), y = 95, label = expression(italic("Acer rubrum")), 
         parse = T, size = 3, colour = "#E41A1C") + #add text to figure
  annotate("text", x = as.Date("2020-3-1"), y = 10, label = expression(italic("Banksia integrifolia")), 
           parse = T, size = 3, colour = "#377EB8") + #add text to figure  
  annotate("text", x = as.Date("2020-3-1"), y = 16, label = expression(italic("Liriodendron tulipifera")), 
           parse = T, size = 3, colour = "#4DAF4A") + #add text to figure
  annotate("text", x = as.Date("2020-1-28"), y = 81, label = expression(italic("Syzygium floribundum")), 
           parse = T, size = 3, colour = "#984EA3") + #add text to figure
  theme_classic() +
  theme(legend.position = "none")

fig3b

#ggsave("Fig3b_R1.pdf", height = 5, width = 5)
# Simulate basic survival curves by life stage
time3 <- seq(1,19,1)
seed <- c(rep(1,8),c(1,0.9,0.7,0.6,0.5,0.45,0.4,0.4,0.4,0.4,0.35))
seedling <- c(rep(1,6),c(0.9,0.8,0.8,0.7,0.6,0.4,0.25,0.2,0.2,0.15,0.1,0.1,0.1))
immature <- c(rep(1,8),c(1,1,0.9,0.9,0.9,0.85,0.85,0.8,0.7,0.65,0.65))
adult <- c(rep(1,8),c(1,1,1,1,1,1,1,1,0.8,0.7,0.7))
demodat_c <- data.frame(cbind(time3, seed, seedling, immature, adult))
demodat2 <- pivot_longer(demodat_c, cols = c(seed, seedling, immature, adult), names_to = "stage")
demodat2$treatment <- "control"

seed <- c(rep(1,7),c(0.9,0.8,0.8,0.6,0.5,0.4,0.3,0.3,0.3,0.3,0.3,0.25)*0.5)
seedling <- c(rep(1,6),c(0.8,0.7,0.7,0.5,0.2,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1)*0.3)
immature <- c(rep(1,8),c(1,0.7,0.6,0.6,0.6,0.55,0.55,0.5,0.45,0.45,0.45)*0.7)
adult <- c(rep(1,10),c(1,0.8,0.8,0.8,0.75,0.7,0.6,0.55,0.55)*0.8)
demodat_h <- data.frame(cbind(time3, seed, seedling, immature, adult))
demodat3 <- pivot_longer(demodat_h, cols = c(seed, seedling, immature, adult), names_to = "stage")
demodat3$treatment <- "heat"
demodat4 <- rbind(demodat2,demodat3)
demodat4$stage <- factor(demodat4$stage, levels = c("seed","seedling","immature","adult"))

fig3d <- 
ggplot(demodat4, aes(x = time3, y = value, colour = stage, lty = treatment)) + 
  geom_line() +
  scale_colour_manual(values = cbp2[c(2,3,4,6)], name = "") +
  scale_linetype_manual(values = c(1,2), labels = c("Tolerant", "Sensitive"), name = "") +
  labs(y = "Survival probability to reach \n next life stage or to reproduce", 
       x = "Cumulative thermal load") +
  theme_classic() + theme(axis.text = element_blank(), legend.position = "bottom") + 
  guides(colour = guide_legend(nrow = 2, byrow = TRUE), linetype = guide_legend(nrow = 2, byrow = TRUE))

fig3d

#ggsave("Fig3d_R1.pdf", width = 5, height = 5)
# Simulate temperature and multi-stressor data
simTemp <- c(30:70)
simTime <- seq(0,2.5,0.061)
simH <- seq(55,35,-0.5) + rnorm(seq(55,35,-0.5),0,0.5)
simA <- seq(50,30,-0.5) + rnorm(seq(50,30,-0.5),0,0.5)
simB <- seq(55,25,-0.74) + rnorm(seq(55,25,-0.74),0,0.5)
simP <- seq(50,20,-0.74) + rnorm(seq(50,20,-0.74),0,0.5)
simAB <- seq(60,28,-0.8) + rnorm(seq(60,28,-0.8),0,0.7)
simTA <- seq(55,29,-0.65) + rnorm(seq(55,29,-0.65),0,0.5)
simTB <- seq(55,25,-0.75) + rnorm(seq(55,25,-0.75),0,0.5)
simTC <- seq(52,12,-1)  + rnorm(seq(52,12,-1),0,1)

simdata <- data.frame(simTime, simTemp, simH, simA, simB, simP, simAB, simTA, simTB, simTC)
simdata <- simdata %>% filter(row_number() %% 2 != 1)

simdataA <- simdata %>% pivot_longer(cols = c(simH, simA, simB, simAB), names_to = "sim_type")
simdataA1 <- simdata %>% pivot_longer(cols = c(simH, simP), names_to = "sim_type")
simdataB <- simdata %>% pivot_longer(cols = c(simH, simTA, simTB, simTC), names_to = "sim_type")

sim_mod <- lm(value ~ simTime * sim_type, simdataA1)
simH_mod <- lm(value ~ simTime, subset(simdataA1, sim_type=="simH"))
simP_mod <- lm(value ~ simTime, subset(simdataA1, sim_type=="simP"))

summary(sim_mod)
## 
## Call:
## lm(formula = value ~ simTime * sim_type, data = simdataA1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9067 -0.2909 -0.1095  0.2975  1.0898 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           55.1223     0.2187  252.06   <2e-16 ***
## simTime               -8.1656     0.1553  -52.59   <2e-16 ***
## sim_typesimP          -5.2454     0.3093  -16.96   <2e-16 ***
## simTime:sim_typesimP  -3.8012     0.2196  -17.31   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4885 on 36 degrees of freedom
## Multiple R-squared:  0.9972, Adjusted R-squared:  0.997 
## F-statistic:  4266 on 3 and 36 DF,  p-value: < 2.2e-16
summary(simH_mod)
## 
## Call:
## lm(formula = value ~ simTime, data = subset(simdataA1, sim_type == 
##     "simH"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89227 -0.23981 -0.02263  0.28588  0.84318 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  55.1223     0.1995  276.34   <2e-16 ***
## simTime      -8.1656     0.1416  -57.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4456 on 18 degrees of freedom
## Multiple R-squared:  0.9946, Adjusted R-squared:  0.9943 
## F-statistic:  3323 on 1 and 18 DF,  p-value: < 2.2e-16
summary(simP_mod)
## 
## Call:
## lm(formula = value ~ simTime, data = subset(simdataA1, sim_type == 
##     "simP"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9067 -0.2947 -0.1378  0.3041  1.0898 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  49.8769     0.2363  211.04   <2e-16 ***
## simTime     -11.9668     0.1678  -71.31   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.528 on 18 degrees of freedom
## Multiple R-squared:  0.9965, Adjusted R-squared:  0.9963 
## F-statistic:  5085 on 1 and 18 DF,  p-value: < 2.2e-16
fig3e <-
ggplot(data = simdataA, aes(x = simTime, y = value, colour = sim_type)) +
  geom_smooth(aes(y = value), method = "lm", formula = y~x, fullrange = T, se = F) +
  labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
  coord_cartesian(ylim = c(41,60)) + 
  scale_x_continuous(expand = c(0, 0), limits = c(0,2.5)) +
  scale_colour_manual(values = c(cbp2[8],cbp2[3],cbp2[4],cbp2[1]),
                      labels = c(bquote(Stressor~altering~italic(CT[max1h])),
                                 bquote(Stressor~altering~italic(CT[max1h])~and~italic(z)),
                                 bquote(Stressor~altering~italic(z)),
                                "Heat stress only"), name = "") +
  theme_classic() + theme(legend.position = "bottom") + 
  guides(colour = guide_legend(nrow = 2, byrow = TRUE), linetype = guide_legend(nrow = 2, byrow = TRUE))

fig3e

#ggsave("Fig3e_R1.pdf", height = 5, width = 5)


fig3f <-
ggplot(data = simdataB, aes(x = simTime, y = value, colour = sim_type)) +
  geom_smooth(aes(y = value, linetype = sim_type), method = "lm", formula = y~x, fullrange = T, se = F) +
  labs(y = bquote(Temperature~(degree*C)), x = bquote(log[10]~Time~(h))) +
  coord_cartesian(ylim = c(28,60)) + scale_x_continuous(expand = c(0, 0), limits = c(0,2.5)) +
  scale_colour_manual(values = c(cbp2[1],cbp2[2],cbp2[7],cbp2[7]),
                      labels = c("Heat stress only",
                                 "Heat stress + Stressor A",
                                 "Heat stress + Stressor B",
                                 "Heat stress + Stressors A + B"), name = "") +
  scale_linetype_manual(values = c(1,3,3,1),
                      labels = c("Heat stress only",
                                 "Heat stress + Stressor A",
                                 "Heat stress + Stressor B",
                                 "Heat stress + Stressors A + B"), name = "") +
  theme_classic() + theme(legend.position = "bottom") + 
  guides(colour = guide_legend(nrow = 2, byrow = TRUE), linetype = guide_legend(nrow = 2, byrow = TRUE))

fig3f

#ggsave("Fig3f_R1.pdf", height = 5, width = 5)
# Generate tolerance landscape for parameterising models
# Drosophila suzukii

data <- read.csv("Dsuzukii_data_Orsted2024.csv")
head(data)
##           id temp     time   sex rep prod dead t_coma
## 1 38-300-F-1   38 48.30624 FALSE   1    0    1  14.50
## 2 38-300-F-2   38 48.30624 FALSE   2    0    1  20.80
## 3 38-300-F-3   38 48.30624 FALSE   3    0    1  13.42
## 4 38-300-F-4   38 48.30624 FALSE   4    0    1  16.25
## 5 38-300-F-5   38 48.30624 FALSE   5    0    1  12.50
## 6 38-300-F-6   38 48.30624 FALSE   6    0    1   7.50
tl <- tolerance.landscape(data$temp, data$time)

Tb <- data$temp
time <- data$time
surv.time <- tl$S[ ,2]
surv.pct <- seq(100, 0, length.out = length(c(0, surv.time)))
par(mfrow = c(1, 1), mar = c(4.5, 4, 1, 1))
plot(c(0, surv.time), surv.pct, type = 'l', ylab = 'survival, %', xlab = 'time')

Tb.mn <- tl$ta.mn
z <- 3.27 # tl$z
slope <- -1/z
intercept <- 11.092 # -tl$ctmax * slope for hours. # Set to 11.092 to replicate ctmax.4h value from Ørsted et al 2024
plot(seq(30, 41), 10 ^ (seq(30, 41) * slope + intercept), 
     ylab = 'survial time, hours', xlab = 'body temperature, C', pch = 16, type = 'b')

ctmax.1h <- (log10(1) - intercept) / slope # temperature causing 50% mortality in 1 hour
ctmax.4h <- (log10(4) - intercept) / slope
ctmax.24h <- (log10(24) - intercept) / slope # temperature causing 50% mortality in 1 day
R0 <- 1
k <- log(10) / z


# Set Arrhenius model and repair rate parameters        
TA <- 14065 # Arrhenius temperature, K
TL <- 10.5 + 273.15 # Arrhenius temperature lower threshold, K
TH <- 28.5 + 273.15 # Arrhenius temperature upper threshold, K
TAL <- 50000 # Arrhenius temperature lower, K
TAH <- 100000 # Arrhenius temperature upper, K
TREF <- 20 + 273.15 # reference temperature, K
Tbs <- seq(0, 50) # sequence of temperatures over which to model, deg C
kdot_ref <- 0.00 # repair rate % per minute at ref temp 20 deg C

# Plot the thermal response curve for repair
repair1 <- ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot_ref) 
par(mfrow = c(1, 1), mar = c(4.5, 4, 1, 1))
plot(Tbs, repair1, type = 'l', ylab = 'repair, % per min', xlab = 'body temperature')
abline(v = 28.5, col = 2)

#Tb <- seq(30, 37,1/(60 * 24))
#Tb <- rep(37, 60 * 24)
Tb <- rep(ctmax.1h, 60 * 24)
shift <- 10 ^ ((Tb.mn - Tb) / z) 
# minutes to shift survival curve by at each minute given temperature at that minute    
time.rel <- 0
result <- matrix(data = 0, nrow = length(Tb) - 1, ncol = 2)
repair <- ArrFunc5(Tb, TA, TAL, TAH, TL, TH, TREF, kdot_ref) 

for(i in 1:(length(Tb) - 1)){
#for(i in 130:150 - 1){
  # adjust time for the survival probability vs time relationship to account 
  # for current temperature and then get % alive after the current time increment
  alive <- try(approx(c(0, shift[i] * surv.time), surv.pct, 
                      xout = time.rel + 1, ties = "ordered")$y, silent = TRUE)
  alive <- min(100, alive + repair[i]) # add repair, but still catch NA
  if(is.na(alive)){ # cap at 100% surviving
    alive <- 100
    break
  }
  result[i, 1]<- alive # save current survival fraction
  # now interpolate to get the time that goes with the survival %
  time.rel <- try(approx(surv.pct, c(0, shift[i + 1] * surv.time),
                         xout = alive)$y, silent = TRUE)
  result[i, 2] <- time.rel # save time since start
}

plot(result[, 2], result[, 1], type = 'l', ylab = "% surviving", xlab = "time, min")

# Integrate damage (experimental end point of death/coma gives lethal dose)
L <- matrix(0, nrow = length(Tb))

for(i in 2:(length(Tb))){
  if(Tb[i] > ctmax.24h){
    dLdt <- exp(k * (Tb[i] - ctmax.24h) - 1)
  }else{
    dLdt <- 0
  }
  L[i] <- L[i-1] + dLdt
}

Ld <- L[which(result[, 1]==0)[1]]
Ld
## [1] 1748.163

Spatial distribution simulation

# day of year to start simulation
DOY <- 1
sim.length <- 7 # n days to simulate

# read in air temps and filter for month, ctmax
tair.120cm <- brick("TA120cm_2017.nc")
## Error in charToDate(x) : 
##   character string is not in a standard unambiguous format
tair.120cm.sub <- tair.120cm[[(DOY * 24 + 1):(DOY * 24 + sim.length * 24)]]
tair.ctmax <- tair.120cm.sub / 10
tair.ctmax[tair.ctmax < ctmax.1h] <- 0
tair.ctmax[tair.ctmax > ctmax.1h] <- 1
sum.sub <- calc(tair.ctmax, function(x){sum(x)})
max.sub <- calc(tair.120cm.sub, function(x){max(x)})

ct.range <- sum.sub
ct.range[ct.range == 0] <- -1
ct.range[ct.range >= 0] <- 0
ct.range[ct.range == -1] <- 1

# create grid of points
nc <- nc_open("TA120cm_2017_time.nc")
lon <- ncvar_get(nc, "longitude") # extract longitude values
lat <- ncvar_get(nc, "latitude") # extract latitude values
sites <- expand.grid(lon, lat) # 3618 sites to iterate over each day

# Loop across sites and calculate damage (first with no repair) 
# for dynamic tolerance landscapes (Rezende approach) and dynamic ctmax (Jørgensen approach) using the 
# Drosophila suzukii female fecundity parameterised data plotted over Australia.
# Note: Simulation may take well over 1 hour to run for each level of repair (0 or 0.003)

kdot_ref <- 0 ### change to 0.02 to simulate repair, or 0 for no repair ###
all.results <- matrix(nrow = nrow(sites), ncol = 4)

par(mfrow = c(1, 1), mar = c(4.5, 4, 1, 1))
plot(ct.range, zlim = c(0.5, 1), col = 'grey', main = "range based on 120 cm air temperature, deg C", 
     ylab = "latitude", xlab = "longitude")
ozmap(x = "country", add = TRUE)

Code for running spatial distribution simulation (not run in markdown)

# parameters referring to survival or surv can be any sublethal measure, this is a generic function term

kdot_ref <- 0

for(ii in 1:nrow(sites)){ # start site loop
  loc <- c(sites[ii, 1], sites[ii, 2])
  tair.120cm_loc <- extract.nc("TA120cm_2017_time.nc", loc)[(DOY * 24 + 1):(DOY * 24 + sim.length * 24), ]
  if(!is.na(tair.120cm_loc$data[1])){ # check if land or sea
    if(is.na(all.results[ii, 1])){
      variable.temp <- tair.120cm_loc$data / 10
      if(max(variable.temp) < tl$ctmax){ # don't do sites that won't survive one minute
        time.min <- spline(variable.temp, n = length(variable.temp - 1) * 60)$x - 1 
        Tb.min <- spline(variable.temp, n = length(variable.temp - 1) * 60)$y
        
        # Daily survival probability
        # dynamic.landscape function (function returns NA for survival probability < 0) 
        # Loop every 1440 min (= 24 h * 60 min)
        day <- rep(1:sim.length, each = 1440)
        surv.out <- seq(100, 0, length.out = length(c(0, surv.time))) # vector to hold survival
        # plot(surv.time, surv.pct[-1], ylab = 'surviving %', xlab = 'time, min', type = 'l')
        final <- matrix(data = 100, nrow = sim.length, ncol = 1)
        final2 <- final
        alive <- 100
        last.alive <- alive
        damage <- 0
        result2 <- matrix(data = 0, nrow = length(tair.120cm_loc$data) * 60, ncol = 2)
        for(xx in 1:sim.length){ # loop through days
          Tb <- Tb.min[xx == day]
          shift <- 10 ^ ((Tb.mn - Tb) / z) # minutes to shift survival curve by at each minute  
          time.rel <- 0
          result <- matrix(data = 0, nrow = length(Tb) - 1, ncol = 2)
          repair <- ArrFunc5(Tb, TA, TAL, TAH, TL, TH, TREF, kdot_ref)
          #plot(seq(0, 1440-1), Tb, ylab = 'body temperature', xlab = 'time, min', type = 'l')
          #plot(seq(0, 1440-1), repair, ylab = '% repair / min', xlab = 'time, min', type = 'l')
          #plot(seq(0, 1440-1), shift, ylab = 'survival shift, min', xlab = 'time, min', type = 'l')
          # count <- 0
          if(xx == 1){
            starti = 2
            result[1, 1] <- 100
          }else{
            starti = 1
          }
          for(i in starti:(length(Tb) - 1)){
            alive <- try(approx(c(0, shift[i] * surv.time),
                                surv.pct, xout = time.rel + 1, ties = "ordered")$y,
                         silent = TRUE)
            if(is.na(alive)){
              if(last.alive != 100){
               alive <- 0
              }else{
               alive <- 100
              }
            }
            if(alive == 0){
             break
            }
            alive <- min(100, alive + repair[i]) # add repair, but still catch NA
            # modify kdot_ref based on alive. If alive > 99, kdot_ref stays the same. 
            # If alive < 99, kdot_ref reduces as a function of alive
            kdot_ref1 <- ifelse(alive < 99, kdot_ref*(kdot_ref^(3/(alive))), kdot_ref)
            alive <- alive + ArrFunc5(Tb[i], TA, TAL, TAH, TL, TH, TREF, kdot_ref1)
            alive <- ifelse(alive > 100, 100, alive)
            last.alive <- alive
            result2[(xx - 1) * 24 * 60 + i, ] <- c(i, alive)

            result[i, 1] <- alive
            time.rel <- try(approx(surv.pct, c(0, shift[i + 1] * surv.time), xout = alive)$y, silent = TRUE)
            result[i, 2] <- time.rel
          }
          
          cat(xx, '\n')
          final[xx] <- min(result[, 1], na.rm = TRUE)
          
          # Jørgensen et al. approach - adapted from equation 2: https://doi.org/10.1038/s41598-021-92004-6
          L <- matrix(Ld, nrow = length(Tb))
          L[1] <- damage
          for(i in 2:(length(Tb))){
            if(Tb[i] > ctmax.24h){
              dLdt <- exp(k * (Tb[i] - ctmax.24h) - 1)
            }else{
              dLdt <- 0
            }
            L[i] <- max(0, L[i-1] + dLdt - repair[i] / 100 * Ld)
            if(L[i] > Ld){ # dead
              L[i] <- Ld
              break
            }
          }
          #points(L/Ld * 100, col = 'purple', type = 'l')
          # work out damage end point for carry over to next day
          damage.end <- tail(L, 1)
          if(damage.end > Ld){
           damage <- Ld
          }
          final2[xx] <- max(L, na.rm = TRUE)
        } # end of day loop
        #plot(100 - result2[, 2], type = 'l')
        #plot(100 - final, type = 'h')
        #points(final2 / Ld * 100, type = 'p', pch = 16, col = 2, ylim = c(0, 100))
        #plot(100-final, final2 / Ld * 100, pch = 16, ylim = c(0, 100), xlim = c(0, 100), ylab = 'Jorgensen', xlab = 'Rezende')
        #abline(0, 1)
        points(loc[1], loc[2], cex = sum(final[,1])/30/100, pch = 16, col = 'blue')
        if(min(final[,1]) / 75 == 0){
          points(loc[1], loc[2], pch = 4)
        }else{
          points(loc[1], loc[2], cex = min(final[,1]) / 75, pch = 16, col = 'blue')
        }
        # write.table(c(loc, min(final[,1])), file = "Rezende.out.rev0.csv", sep = ",",
        #            col.names = F, qmethod = "double", append = T)
        # write.table(c(loc, max(final2[,1])/Ld * 100), file = "Jorgensen.out.rev0.csv", sep = ",",
        #            col.names = F, qmethod = "double", append = T)
        all.results[ii, ] <- c(loc, min(final[, 1]), max(final2[, 1]) / Ld * 100) 
        cat('site', ii, '\n')
      }else{
        points(loc[1], loc[2], pch = 4)
      }
    }
  }
}

all.results2 <- as.data.frame(na.omit(all.results))
colnames(all.results2) <- c('lon', 'lat', 'surv', 'dose')
plot(ct.range, zlim = c(0.5, 1), col = 'grey',  ylab = 'latitude', xlab = 'longitude', legend = FALSE)
points(all.results2[, 1:2], cex = 1 - all.results2$dose / 75, pch = 16, col = 'blue')
plot(ct.range, zlim = c(0.5, 1), col = 'grey',  ylab = 'latitude', xlab = 'longitude', legend = FALSE)
points(all.results2[, 1:2], cex = all.results2$surv / 120, pch = 16, col = 'blue')
all.results.no.repair <- all.results2


#write.csv(all.results.no.repair, 'dsuz.results.no.repair0.csv')
#write.csv(all.results.with.repair, 'dsuz.results.repair0.02.csv')
dsuz.results.no.repair0 <- read.csv("dsuz.results.no.repair.pro0.csv")[, -1]
dsuz.results.repair2 <- read.csv("dsuz.results.repair.pro0.02.csv")[, -1]
rezende.no.repair <- read.csv("Rezende.out.pro0.csv", header = F, col.names = c("Var","Value"))
jorgensen.no.repair <- read.csv("Jorgensen.out.pro0.csv", header = F, col.names = c("Var","Value"))

head(rezende.no.repair)
##   Var     Value
## 1   1 132.11900
## 2   2 -11.40100
## 3   3  99.41495
## 4   1 142.31900
## 5   2 -11.40100
## 6   3  99.43911
head(jorgensen.no.repair)
##   Var   Value
## 1   1 132.119
## 2   2 -11.401
## 3   3   0.000
## 4   1 142.319
## 5   2 -11.401
## 6   3   0.000
rezende.no.repair_wide <- pivot_wider(rezende.no.repair, names_from = Var, values_from = Value)
rezende.no.repair_wide <- cbind(data.frame(rezende.no.repair_wide$`1`, rezende.no.repair_wide$`2`, rezende.no.repair_wide$`3`))
names(rezende.no.repair_wide) <- c("lat", "lon", "r_surv")
head(rezende.no.repair_wide)
##       lat     lon   r_surv
## 1 132.119 -11.401 99.41495
## 2 142.319 -11.401 99.43911
## 3 132.719 -12.001 77.52720
## 4 133.319 -12.001 49.95959
## 5 133.919 -12.001 66.39607
## 6 142.319 -12.001 88.41603
jorgensen.no.repair_wide <- pivot_wider(jorgensen.no.repair, names_from = Var, values_from = Value)
jorgensen.no.repair_wide <- cbind(data.frame(jorgensen.no.repair_wide$`1`, jorgensen.no.repair_wide$`2`, jorgensen.no.repair_wide$`3`))
names(jorgensen.no.repair_wide) <- c("lat", "lon", "j_surv")
head(jorgensen.no.repair_wide)
##       lat     lon    j_surv
## 1 132.119 -11.401  0.000000
## 2 142.319 -11.401  0.000000
## 3 132.719 -12.001 16.770258
## 4 133.319 -12.001 42.385432
## 5 133.919 -12.001 27.874079
## 6 142.319 -12.001  8.570505
no.repair.dat <- as.data.frame(cbind(rezende.no.repair_wide, (100-jorgensen.no.repair_wide[,3])))
names(no.repair.dat) <- c("lat", "lon", "r_surv", "j_surv")
no.repair.dat$diff <- no.repair.dat$r_surv - no.repair.dat$j_surv
no.repair.dat$diff.pos <- ifelse(no.repair.dat$diff < 0, NA, no.repair.dat$diff)
no.repair.dat$diff.neg <- ifelse(no.repair.dat$diff > 0, NA, no.repair.dat$diff)

tair.brks <- seq(20, 50, 5) # breaks for the colour scheme
tair.col <- colorRampPalette(c("blue", "yellow", "red"))(length(tair.brks) - 1) # colour
# Multi-panel plot
par(mfrow = c(2,2), mar = c(3,1,3,1), oma = c(1,1,3,1))

# Maximum air temp over the simulation
plot(max.sub / 10, breaks = tair.brks, col = tair.col, ylab = 'latitude', xlab = 'longitude', axes = F, box = F,
     main = bquote(atop(bold(a)~"Maximum air temperature"~(120~cm~","~degree*C),
                        Productivity~probability~("%")~italic(CT["max1h"])~model)), alpha = 0.3)
plot(max.sub / 10, zlim = c(20, ctmax.1h), breaks = tair.brks, col = tair.col, 
     ylab = 'latitude', xlab = 'longitude', add = T, legend = F)
ozmap(x = "country", add = TRUE)

# Rezende model
plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
     main = bquote(atop(bold(b)~Productivity~probability~("%"),~without~repair~model)))
points(dsuz.results.no.repair0[, 1:2], cex = dsuz.results.no.repair0$surv/120, pch = 16, col = cbp2[6])
ozmap(x = "country", add = TRUE)
legend(x = 124, y = -32, legend = paste(seq(20, 100, 20), "%"), pt.cex = seq(20, 100, 20)/75, pch = 16, col = cbp2[6], bty = 'n')

# Jørgensen model

# plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
#      main = bquote(atop(bold(b)~Survival~probability~("%")~without~repair,Jørgensen~italic(et~al.)~dynamic~italic(CT[max])~model)))
# points(no.repair.dat[, 1:2], cex = no.repair.dat$j_surv/120, pch = 16, col = cbp2[3])
# ozmap(x = "country", add = TRUE)
# legend(x = 126, y = -32, legend = paste(seq(20, 100, 20), "%"), pt.cex = seq(20, 100, 20)/75, pch = 16, col = cbp2[3], bty = 'n')


# Rezende model with repair at kdot=0.02
plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
     main = bquote(atop(bold(c)~Productivity~probability~("%"),~with~high~repair~rate~model)))
points(dsuz.results.repair2[, 1:2], cex = dsuz.results.repair2$surv / 120, pch = 16, col = cbp2[4])
ozmap(x = "country", add = TRUE)
legend(x = 124, y = -32, legend = paste(seq(20, 100, 20), "%"), pt.cex = seq(20, 100, 20)/75, pch = 16, col = cbp2[4], bty = 'n')

# Diff between no repair and repair at 0.02
dsuz.results.repair2$diff <- dsuz.results.repair2$surv - dsuz.results.no.repair0$surv
range(dsuz.results.repair2$diff)
## [1] -10.59975  10.67527
plot(ct.range, zlim = c(0.5, 1), col = 'grey', ylab = 'latitude', xlab = 'longitude', legend = FALSE, axes = F, box = F,
     main = bquote(atop(bold(d)~Productivity~probability~difference~("%"),~with~repair~"-"~without~repair~models)))
points(dsuz.results.repair2[, 1:2], cex = dsuz.results.repair2$diff/15, pch = 16, col = cbp2[8])
ozmap(x = "country", add = TRUE)
legend(x = 124, y = -32, legend = paste(seq(5, 25, 5), "%"), pt.cex = seq(5, 25, 3)/15, pch = 16, col = cbp2[8], bty = 'n')

# Save figure as approximately 8 x 8 inch